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Abstract. This article is concerned with numerical methods to approximate effective 
coefficients in stochastic homogenization of discrete linear elliptic equations, and their 
numerical analysis — which has been made possible by recent contributions on quantita- 
tive stochastic homogenization theory by two of us and by Otto. This article makes the 
connection between our theoretical results and computations. We give a complete pic- 
ture of the numerical methods found in the literature, compare them in terms of known 
(or expected) convergence rates, and study them numerically. Two types of methods 
are presented: methods based on the corrector equation, and methods based on random 
walks in random environments. The numerical study confirms the sharpness of the anal- 
ysis (which it completes by making precise the prefactors, next to the convergence rates), 
supports some of our conjectures, and calls for new theoretical developments. 
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1. Introduction 

There have always been strong connections between stochastic homogenization of linear 
elliptic PDEs and random walks in random environments (RWRE). In the first case, 
the central object is the random elliptic operator which can be replaced on large scales 
by a deterministic elliptic operator characterized by the so-called homogenized matrix. 
In the second case, the distribution of the rescaled random walk approaches that of a 
Brownian motion, whose covariance matrix is deterministic (and corresponds to twice the 
homogenized matrix). 

Both approaches allow one to devise numerical methods to approximate the homoge- 
nized coefficients. Until very recently, stochastic homogenization in general was only a 
qualitative theory, so that there was no quantitative analysis of these numerical methods. 
In a series of papers by two of us, Otto, and Neukamm, we developed a quantitative theory 
of stochastic homogenization |GQ1H IGQ12al IMolU IGNOaj which allowed us to make a 
quantitative numerical analysis of these methods |G112| IGMl2al IGM12b] IGNObj . 

In the case of stochastic homogenization of discrete elliptic equations with stationary 
random conductances of finite correlation length, and the corresponding random walk on 
Z d , it is now possible to give a complete quantitative picture of numerical methods to 
approximate the homogenized coefficients. 
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The aim of this article is twofold. First we address both approaches of stochastic 
homogenization together: the corrector equation for the PDE approach, and the RWRE for 
the probabilistic approach. Since the mathematics behind both approaches are somewhat 
different, they are rarely presented together, whereas for the quantitative analysis both 
can be needed (see for instance |GM12al IGM12bj ). A crucial role is played by spectral 
analysis to connect both worlds, and for the quantitative analysis. We shall present and 
compare several ways of approximating homogenized coefficients, based on the corrector 
equation and on the RWRE. For each strategy we provide the current state of the art 
of the numerical analysis, and motivate conjectures when the analysis is not complete. 
Second, we report a systematic study of these numerical strategies. Our numerical results 
confirm the sharpness of the numerical analysis, support some of our conjectures, and call 
for new theoretical insight. 

For completeness, and in order to make this article accessible to probabilists, analysts, 
and numerical analysts, we provide a short comprehensive review of the well-known qual- 
itative results of stochastic homogenization, as well as informal arguments for the quan- 
titative results. We believe that this article, which sets and compares on rigorous ground 
various numerical methods on the prototypical example of stochastic homogenization of 
discrete elliptic equations, will also be valuable to practitioners. 

2. Stochastic homogenization: corrector equation and RWRE 

We start with the description of the discrete set of diffusion coefficients, first present the 
discrete elliptic point of view, and then turn to the random walk in random environment 
viewpoint. The aim of this section is to introduce a formalism, and give an intuition on 
both points of view. 

The results recalled here are essentially due to Papanicolaou and Varadhan |P V79j . 
Kozlov }Ko87j , and Kipnis and Varadhan |KV86j ). 

We present the results in the case of independent and identically distributed (i.i.d.) con- 
ductivities, although everything in this section remains valid provided the conductivities 
lie in a compact set of (0, +oo), are stationary, and ergodic. In particular, we shall apply 
this theory (and its quantitative counterpart) to coefficients which have finite correlation 
length. Yet, for the sake of clarity, we stick to i.i.d. coefficients in this presentation. 

2.1. Random environment. We say that x,y in Z d are neighbors, and write x ~ y, 
whenever \y — x\ = 1. This relation turns Z rf into a graph, whose set of (non-oriented) 
edges is denoted by B. We now define the conductances on B, and their statistics. 

Definition 2.1 (environment). Let < a ^ (3 < +oo, and Q = [a,/3] B . An element 
uj = (o; e ) e gB of £1 is called an environment. With any edge e = (x,y) £ B, we associate 
the conductance w^y) := w e (by construction ^( x ,y) = ^(y.x))- Let v be a "probability 
measure on [a, 0\. We endow f2 with the product probability measure P = v®^. In other 
words, if lo is distributed according to the measure P, then (w e ) eG B ar & independent random 
variables of law v. We denote by L 2 (Q) the set of real square integrable functions on Q, 
for the measure ¥, and write {■) for the expectation associated with P. 

We then introduce the notion of stationarity. 
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Definition 2.2 (stationarity). For all z £ 7L d , we let 9 Z : ft — > SI be such that for all to £ 
and £ B, (6> 2 = W( I + z . y + x ) • TViis defines an additive action group {0 z } zeZ d 

on ft which preserves the measure P, and is ergodic for P. 

We say that a function f : Q, x 7L d — > E is stationary if and only if for all x, z £ Z d and 
P-almost every uj £ Q, 

f(x + z,u) = f(x,9 z u). 
In particular, with all f £ L 2 (S1), one may associate the stationary function (still denoted 
by f) : Z d x E, (x, w) i — ^ /(^ uj). In what follows we will not distinguish between 
f £ L 2 (S1) and its stationary extension on Z d x SI. 

2.2. Corrector equation. We associate with the conductances on B a conductivity ma- 
trix on 7h d . 

Definition 2.3 (conductivity matrix). Let £1, P, and {6 z } Z £Z d be as in Definitions \2.1\ 
and \2.2\ The stationary diffusion matrix A : Z, d x ft — > .A4d(E) is defined by 

A : (x,w) i-> diag(w (a . ja . + e i ) ! ■ ■ ■ > w (s,s+e d ))- 
For each uj £ O, we consider the discrete elliptic operator L defined by 

L = -V*-A(-,w)V, (2.1) 

where V and V* are the forward and backward discrete gradients, acting on functions 

u : Z d ->• E by 



V-u(x) := 



u(x + ei) — u(x) 



u(x + e d ) 



u[x 



, V*u{x) 



u{x) — u{x — ei) 



u{x) — u{x — ed) 



(2.2) 



and we denote by V*- the backward divergence. In particular, for all u 



Lu : 



z i y u - 



(z,z')(u(z) ~ U{z')). 



(2.3) 



The standard stochastic homogenization theory for such discrete elliptic operators (see 
for instance |Ku83| . |Ko87] ) ensures that there exist homogeneous and deterministic co- 
efficients A\ lom such that the solution operator of the continuum differential operator 
—V • j4h om V describes P-almost surely the large scale behavior of the solution operator of 
the discrete differential operator —V* • A{-,uj)V. As for the periodic case, the definition 
of ^4hom involves the so-called correctors. Let £ £ E rf be a fixed direction. The corrector 
cf) : 7L d x ft — > E in the direction £ is the unique solution (in a sense made precise below) 
to the equation 

- V* • A(x,u)(€ + V<f>{x,u)) = 0, x£Z d . (2.4) 
The following lemma gives the existence and uniqueness of this corrector <p. 

Lemma 2.4 (corrector). Letfl, P, {Oz} z eZ d > an d A be as in Definitions \2.1\ POl and \2.?A 
Then, for all £ £ E d , there exists a unique measurable function 4> '■ Z d x S7 — > E such that 
(p(0, •) = 0, S/cp is stationary, (Vcft) = 0, and 4> solves (|2.4p P-almost surely. Moreover, 
the symmetric homogenized matrix ^4hom is characterized by 

Z ■ Aom£ = <(£ + ■ A{£ + V0)> . (2.5) 
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The standard proof of Lemma [231 makes use of the regularization of (j2.4|) by a zero-order 
term \i > 0: 

/j$p(x,u)-V*-A(x,u)(Z + V<f>p(x,u})) = Q, x£Z d . (2.6) 

Lemma 2.5 (regularized corrector). Let Q, F, {9 z } z ^2, d > an d A be as in Definitions {2H\ 
\2.2l and \2.!A Then, for all \x > and £ G IR rf , i/iere exists a unique stationary function 
(f>n G L 2 (TL) with = which solves (|2.6|) IP -almost surely. 

To prove Lemma 12.51 we follow [PV79] . and introduce difference operators on L 2 {Q): 
for all u G L 2 (0), we set 



Du(w) := 



u(6 ei uj) — u{ui) 



u(6 ed Uj) - U(tj) 



u{oj) — u(9- ei oj) 



u{lo) - u(9^ ed 0j) 



(2.7) 



These operators play the same roles as the finite differences V and V* — this time for the 
variable to (in other words, they define a difference calculus on L 2 (f2)). They allow us to 
define the counterpart on L 2 (Q) to the operator L of (|2. 1 j) : 

Definition 2.6. Let n, F, {9 z } zeZ d, and A be as in Definitions \2^ WM and\2M We 
define C : L 2 {£1) -> L 2 (VL) by 

£u(u) = - D* -A(oj) D u{lo) 

= y)uo,z(u(u) - u(9 z to)) 

where D and D* are as in (|2.7|) . 

The fundamental relation between L and £ is the following identity for stationary fields 



u : 



x n 



l: for all z£Z and almost every oj G f2, 
Lu(z,uj) = £u(9 z u). 

In particular, the regularized corrector 4>^ is also the unique weak solution in L 2 (f2) to the 
equation 

- D* -A(a;)(£ + D ^(w)) =0, 
and its existence simply follows from the Riesz representation theorem on L 2 (Q). 

The regularized corrector is an approximation of the corrector <f> in the following 
sense: 

Lemma 2.7. Let VL, F, {9 z } zeZ d, and A be as in Definitions [KJ\ POl and HOI For all 

\i > and £ G let <f> and 6e i/ie corrector and regularized corrector of Lemmas \2.J\ 
and\2.5[ Then, we have 



lim (ID. 



D, 



0. 



From the elementary a priori estimates 

<|V^| 2 > = <|D^| 2 > ^ C, <^» 2 > < 

for some C independent of //, we learn that D<^ is bounded in L 2 (Q,M. d ) uniformly in \i, 
so that up to extraction it converges weakly in L 2 (Q,,F, d ) to some random field $ (which 



CONVERGENCE RATES IN STOCHASTIC HOMOGENIZATION 5 

is a gradient). This allows one to pass to the limit in the weak formulations and obtain 
the existence of a field <& = ($1, . . . , &d) G L 2 (Q,R d ) such that for all tp G L 2 (Vt), 

(D ip ■ A(£ + $)) = 0. (2.8) 

Using the following weak Schwarz commutation rule 

Vj, k G {1, . . . , d}, ((D,- V)^ fc ) = ((D fc V)^> 

one may define ^ : 7L d x — >• R such that is stationary, $ = V</>, and </>(0, w) = for 
almost every uj G fi. By definition this function is not stationary. It is a priori not clear 
(and even wrong in dimension d ^ 2) whether there exists some function ip G L 2 (f2) such 
that D-0 = $ (this is a major difference with the periodic case). 

The uniqueness of $ is a consequence of Lemma 12.71 which follows from the fact that 
~D<pn is a Cauchy sequence in L 2 {Q). To prove this, we shall appeal to spectral theory. 

The operator C of Definition 12.61 is bounded, self-adjoint, and non-negative on L 2 {Q). 
Indeed, for all tp,x G L 2 (0), we have by direct computations 

x2\l/2 / AJ 1-51. /.2\V2 



Hence, £ admits a spectral decomposition in L 2 (Q). For all 5 G L 2 (Q.) we denote by e s the 
projection of the spectral measure of C on g. This defines the following spectral calculus: 
for any bounded continuous function ^ : [0, +00) — > R+, 



<(*(£)<?)£> = / *(A)de s (A). 

Let £ G R d with |£| = 1 be fixed, and define the local drift as D = - D* G L 2 (Q). For 
all /U ^ i/ > we have M = (/i + £) _1 t) and tzV = (f + ^C)" 1 !), by the Cauchy-Schwarz 
inequality, 

(|D0 M -D0„| 2 ) < a' 1 ((</>„ - 0„)£(0 M - <t>v)) 

= a" 1 (^£0^) - 2a" 1 {fy^Cfyv) + a" 1 (cf> v £(f> v ) 

= a" 1 + C)- 1 ^ + £) _1 i>) - 2a" 1 (t)(/i + Cy 1 ^ + £) _1 5> 

+a~ x + c)- x c{v + £) _1 o> . 

By the spectral formula with functions 

AAA 



*(A) 



(/i + A) 2 ' ( At + A)(i/ + A)' (1/ + A) 



2 ' 



we obtain 

<|D^-D^> « a- ^ -2 ^ ^ - A) + ¥ ± w ) ^(A) 

r ^2 

since < z/ ^ /i. Since the upper bound is independent of z/, we have proved the claim if 
we can show that it tends to zero as fi vanishes. This is a consequence of the Lebesgue 
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dominated convergence theorem provided we show that 

1 



de (A) < oo. (2.10) 

A 

On the one hand, by the a priori estimate of D 

(^C0 M ) < /3<|D^| 2 > < PC. 
On the other hand, by the same type of spectral calculus as above, we have 

(0 M £0 M ) = f j de a (A). 

Estimate (|2.1U|) then follows from the monotone convergence theorem. This concludes the 
proof of Lemma 12.71 

2.3. Random walk in random environment. We now turn our attention to the prob- 
abilistic aspects of homogenization. This presentation is informal. It aims at being ac- 
cessible to non-specialists of probability theory, and at highlighting the inner similarities 
with the corrector approach of subsection 12.21 



2.3.1. The continuous-time random walk. Let the environment cj be fixed for a while (that 
is, we have picked a realization of the conductivities oo e E [a, 0\, e G B). The random walk 
we wish to define, that we will denote by (X t )t£R + , is a random process whose behavior 
is influenced by the environment. 

To the specialist, it can be defined by saying that it is the Markov process whose 
transition rates are the (u; e ) eg B- The Markov property means that given any time t ^ 0, 
the behavior of the process after time t depends on its past only through its location at 
time t. In other words, the process "starts afresh" at time t given its current location. In 
order to give a complete description of the process, it thus suffices to describe its behavior 
over a time interval [0,t], for some t > 0, no matter how small. As t tends to 0, this 
behavior is given by 



P? [X t = z>] 



tuj Z)Z ' + o(t) if z' ~ z, 

l-E^ z ^,y + o(i) Hz' = z, (2-11) 

o(t) otherwise, 



where P^ is the probability measure corresponding to the walk started at z, that is, 
P^[Xq = z] = 1. Equation ()2.1ip is what is meant when it is said that lo Z)Z i is the jump 
rate from z to z' . 

A more constructive way to represent the random walk is as follows. Let the walk be 
at some site z S Z rf at time t, and start an "alarm clock" that rings after a random time 
T which follows an exponential distribution of parameter 



p w («) := J^,*/. (2.12) 
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This means that for any s ^ 0, the probability that T > s is equal to e p "( 2 ) s . When the 
clock rings, the walk chooses to move to one (out of 2d) neighboring site z' with probability 

p(z - z>) := ^f- (2.13) 

and this choice is made independently of the value of T. 

Note that by the Markov property, the fact that the walk has not moved during the 
time interval [t,t + s] should give no information on the time of the next jump. Only 
exponential distributions have this memoriless property. 

Let us see why thus defined, the random walk satisfies (|2.1ip . The probability that the 
clock rings during the time interval [0, t] is 

Since Pu(z) is bounded by 2d/3 uniformly over z, the probability that the walk makes two 
or more jumps is o(t). The probability that it ends up at z' ~ z at time t is thus 

{Pu(z)t + o(t)) p(z z') - o(t) = u Z)Z d + o(t), 

and the probability that it stays still is indeed as in ([2. lip . 

The link between the random walk and the elliptic operators of the previous subsection 
is as follows. We let (Pt)teR+ be the semi-group associated with the random walk, that is, 
for any t ^ and any bounded function / : Z, d — > R, we let 

p t f(z) = Klf&t)], 

where denotes the expectation associated to the probability measure P£, under which 
the random walk starts at z £ 1i d in the environment uj. This is a semi-group since X has 
the Markov property. As we now show, the infinitesimal generator of this semi-group is the 
elliptic operator — L (where L is defined in (|2.ip ). Recall that the infinitesimal generator 
of a semi- group, applied to /, is given by 

lim ( — — — 

t->-0 \ t 

for any / for which the limit exists. From the description (|2.1ip . this limit is easily 
computed. Indeed, 

P t f(z) = E-[/(X t )] 

= {l-p UJ {z)t)f{z)+tY J ^z,z^{z') + o{t) (2.14) 

= f{ z ) + tY J ^{z,z'){f{z')- f(z)) + o(t), (2.15) 



so that by 



lim [ PtfiZ \ = E^')(/(^)-/(^)) = -Lf(z), (2.16) 



and we have identified the infinitesimal generator to be — L, as announced. 

An important feature of this random walk is that the jump rates are symmetric: the 
probability to go from z to z' in an "infinitesimal" amount of time is equal to the probability 
to go from z' to z, as can be seen on (|2.11|) . This may be rephrased as saying that the 
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counting measure on Z (which puts mass 1 to every site) is reversible for the random 
walk. 

If we were running the random walk on a finite graph instead of Z d , one could normalize 
the counting measure to make it into a probability measure, and then interpret reversibility 
as follows. Start the random walk at a point chosen uniformly at random on the graph, 
and let it run up to time t. Reversibility (of the uniform probability measure) is equivalent 
to the fact that the distribution of this trajectory is the same as the distribution of the 
time-reversed path, running from time t to time 0. In particular, since the starting point is 
chosen uniformly at random, the distribution of the walk at any time must be the uniform 
probability measure. Of course, the problem is that on Z d , it does not make sense to say 
that we choose the starting point "uniformly at random", or in other words, we cannot 
normalize the counting measure into a probabability measure. 

In the previous section, we moved from the operator L to its "environmental" version, 
the operator C. This takes a very enlightening probabilistic meaning. Instead of consider- 
ing the random walk itself, we may consider the environment viewed by the particle, which 
is the random process defined as 

t ^ uj(t) := B x{ t)UJ, 

where {0x) x <EZ d are the translations introduced in Definition 12.21 One can convince oneself 
that (u)(t))teR + is a Markov process, and the important point is that its infinitesimal 
generator is precisely —C. A little computation (see for instance [Mold} Proposition 3.1]) 
shows that the reversibility of the counting measure for the random walk translates into 
the reversibility of the measure P for the process (uj(t))t£R + . This reversibility implies 
that (u)(t))t£R + is stationary if we take the initial environment according to P and then 
let the process run, that is, under the measure 

P := PPq, (2.17) 

the so-called annealed measure. Stationarity means that for any positive integer k and 
any s±, . . . , s/- ^ 0, the distribution of the vector (u(si + t), . . . , ui(sk + t)) does not depend 
on t ^ 0. In particular, under the measure Po, the distribution of co(t) is the measure P 
for any t ^ 0. 

To conclude, we argue that from the reversibility follows the fact that the operator 
L is self-adjoint in L 2 (Q) — which we have already seen more directlty in the previous 
subsection. Indeed, the invariance under time reversal implies that, under Po and for 
any t ^ 0, the vectors (w(0), uj(t)) and (u(t),uj(0)) have the same distribution. For any 
bounded functions /, g : fi — > R, we thus have 

E [/(w(0)) g(u(t))] = E [f(co(t)) g(u(0))] , (2.18) 

where we write Eo for the expectation associated to Po- Letting (Vt)teR + be the semi-group 
associated to — C, we arrive at 

E [/M0)) g(u(t))] = (Eg [/(w(0)) g(uj(t))]) 

= (/(a,) E£ \g(u(t))]) = (f V t g) , 

and ([USD thus reads 

(/ Vtg) = (Vtf g) ■ 

This shows that Vt is self-adjoint for any t ^ 0. Passing to the limit as in (|2.16|) . we obtain 
the self-adjointness of C itself. In fact, reversibility of a Markov process and self-adjointness 
of its infinitesimal generator are two sides of the same coin. 



CONVERGENCE RATES IN STOCHASTIC HOMOGENIZATION 



9 



2.3.2. Central limit theorem for the random walk. The aim of this paragraph is to sketch 
a probabilistic argument justifying the following result. 

Theorem 2.8 ( |KV8 6] ) . Under the measure Po and as e tends to 0, the rescaled random 
walk := (y/eX t/e ) te R + converges in distribution (for the Skorokhod topology) to a 

Brownian motion with covariance matrix 2A\ lom , where ^4hom is as in (|2.5p . In other 
words, for any bounded continuous functional F on the space of cadlag functions, one has 



En 



F(X^) —>E[F(B)}, (2-19) 
u 



where B is a Brownian motion started at the origin and with covariance matrix 2j4 hom , 
and E denotes averaging over B. Moreover, for any £ E M d , one has 



H-XtY >2£-A hom £. (2.20) 

t— S- + 00 



Note that the convergence of the rescaled square displacement in (|2.20p does not follow 
from (|2.19p . since the square function is not bounded. 

From now on, we focus on the one-dimensional projections of X. As in subsection 12. 2[ 
we let £ be a fixed vector of M. d . The idea is to decompose £ • X% as 

i ■ X t = M t + Ru (2.21) 

where (Mt)t^o is a martingale and Rt is a (hopefully small) remainder. 
We recall that (Mt)t^o is a martingale under Eq if for all t ^ and s 0, 

E%[M t+s \T t ] = M t , (2.22) 

where Ft is the cr-algebra generated by {M t ,t G [0, t]}. (If M was a gambler's money, one 
should say that he is playing a fair game when (|2.22p holds, since knowing the history up 
to time t, his expected amount of money at time t + s is the amount he has at time t). 

We look for a martingale of the form Mt = x u {Xt) for some function \ u ■ By the Markov 
property of X, (Mt)t^o of this form is a martingale if and only if, for any z £ 7L d and any 
t > 0, 

KbfiXt)] = X"(z), (2.23) 

From (|2.1ip . we learn that 

K\x"{X t )] = ^(z) - sL x "(z) + 0(s 2 ). 

Hence, a necessary condition is that Lx^iz) = 0, and in fact, this condition is also suffi- 
cient. Keeping in mind that we also want the remainder term to be small, we would like 
to be a perturbation of the z i— > £ ■ z, so that a right choice for x^ should be 

where <j) is the corrector of Lemma 12.41 (compare equation (|2.4p to Lx u {z) = 0). The link 
between the corrector equation and the RWRE appears precisely there. We thus define 

M t =t-X t + <t>(X t ,oj) and Rt = -</>(X u u). 

We now argue that the martingale M has stationary increments under the measure Po- For 
simplicity, let us just see that the distribution of M t+S — M t under Po does not depend on s 
(instead of considering vectors of increments). The argument relies on the stationarity of 
the process of the environment viewed by the particle seen above. Given this stationarity, 
it suffices to see that Mt +S — Mt can be written as a function of (ca(t),u)(t + s)) only. This 
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is possible because the increment M( +s — Mf depends only on V<^>, which is a function of 
u) only (contrary to <p itself which is a priori a function of x and u>). 

Martingales are interesting for our purpose since they are "Brownian motions in dis- 
guise". To make this idea more precise, let us point out that any one-dimensional con- 
tinuous martingale can be represented as a time-change of Brownian motion (this is the 
Dubins-Schwarz theorem, see [RY[ Theorem V.1.6]). If (-B*)t>o is a Brownian motion, a 
time-change of it is for instance Mi = B t i . Note that in this example, the time-change 
can be recovered by computing E[Mf] = t 7 . Here, the martingale we consider has jumps 
(since X itself has), which complicates the matter a little, but let us keep this under the 
rug. Intuitively, in order to justify the convergence to a Brownian motion, we should show 
that the underlying time-change grows linearly at infinity. One can check that two incre- 
ments of a martingale over disjoint time intervals are always orthogonal in 1? (provided 
integration is possible). In our case, since Mt has stationary increments, it thus follows 
that letting ct 2 (£) = E [M^], we have 

E [M 2 ] =a 2 (0 t, (2.24) 

so we are in good shape (i.e. on a heuristic level, it indicates that the underlying time- 
change indeed grows linearly). Letting / : z i->- (x ■ £ + 4>(x,io) ■ £) 2 and using (|2.1ip . one 
can write 

E£[(M t -e) 2 ] = Kif(Xt)\ 
= /(0) + t U{z,z>) ■ £ + u) • 2 + o(t) 

= 2t(£ + V0(O, u) ■ • A(0, w)(e + V0(O, u) • +o{t), 
v . ' 

=:t v(ui) 

since 0(0, u) = 0. From the definition of A^ om in (|2.5|) . we thus get 

a 2 ^) = (v) = 2^A hom C. (2.25) 

Provided we can show that the remainder is negligible, this already justifies (|2.20p . In order 
to prove that (y/eM e -i t )t^o converges to a Brownian motion of variance <r 2 (£) as e tends 
to 0, knowing (|2.24p is however not sufficient: one does not recover all the information 
about the time-change by computing Eq[M 2 ] alone. This can be understood from the fact 
that in the Dubins-Schwarz theorem, the time-change that appears can be random itself, 
and is in fact the quadratic variation of the martingale. We will not go into explaining 
what the quadratic variation of a martingale is in general, but simply state that in our 
case, it is given by 

V t = [ v(uj(s)) ds, 
Jo 

and what we should prove is thus that 

t~ l V t = t- 1 [ v(u(s)) ds a\0 = 2£ • A hom £. (2.26) 

One can show that the process {u{t))t^o is ergodic for the measure Po (see for instance 
[MolOl Proposition 3.1]), and thus the convergence in (|2.26p is a consequence of the ergodic 
theorem. Note in passing this surprising fact that the proof of a central limit theorem was 
finally reduced to a law of large numbers type of statement. 
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In order to obtain a central limit theorem for £ X itself (instead of the martingale part), 
we need to argue that the remainder term is small. In order to do so, and following [KV86j . 
we will rely on spectral theory. We recall that the operator C of Definition 12.61 being a 
bounded non-negative self-adjoint operator on L 2 (J7), it admits a spectral decomposition in 
L 2 (Q). Moreover, for any g E L 2 (Q) the characterizing property of the spectral measure e g 
previously defined is that for any continuous function ^ : [0, +oo) — > M + , one has 



{g 9(£)g) = / *(A) de 9 (A). 

We now argue that 

h[R 2 ] = 2 f l —^de,{\) — -> 0, (2.27) 

where Z> is the local drift in the direction £, that is, 

= -v* • A(0, u)Z = - D* -A£ = ■ z. (2.28) 

We start by showing the equality in (|2.27|) . and will later show that the spectral integral 
tends to as t tends to infinity. Recall that 

R t = -</>(X t) u) = 0(0, u) - cf>(X t ,uj), 

and that 0(0, u))—(j)(x, u>) can be obtained as the limit of 0^(0, uj) — cfr^x, uj). To make them 
easier, we do the computations below as if was a stationary 0„. The argument can be 
made rigorous through spectral analysis as in |MolH Theorem 8.1], or using Lemma 12.71 
We expand the square: 

E [R 2 ] =E o [(0(O,o;)) 2 ]+Eo[(0(X t , W )) 2 ]-2Eo[0(O,a;)0(X t ,a;)]. 

Note that, by the simplifying assumption, 

E o [(0(X t ,u;)) 2 ] = E o [(0(O,0 Xt a;)) 2 ] =E o [(0(O,u;(i))) 2 ]. 

Since (u(t))teR + is stationary under Pq, the last expectation is equal to Eo[(0(O, uj)) 2 ] = 
(0 2 ). Since £0 = t), we have 

<0 2 ) = (C-'d C^d) = (d C-H) = J X' 2 de 9 (A). 

For the cross-product, 

E o [0(O,u/)0(X*,u;)] = (0(0, w) E£ [0(0, uj(t))]) , (2.29) 

and Eq [0(0, ui(t))] is the image by the semi-group associated to — C, that is, e~ ll ~ . Using 
also the fact that £0 = 0, we can rewrite the r.h.s. of (|2.29p as 

(0 e- tC 6) ~ 'r- 1 * ~-*r-ir>\ _ f i-2„-tA 



[C-h e'^C-h) = J \- 2 e- tx de B (A), 



and this justifies the equality (|2.27|) . The fact that the spectral integral in (|2.27p tends to 
zero follows from the dominated convergence theorem and (|2.10p . 

Before concluding this section, we point out several differences between the argument 
as we presented it and the original one from |KV86j . Here, we used the existence of 
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the corrector, borrowed from subsection 12. 2\ to construct the martingale. In [KV86J, the 
martingale is constructed directly, by considering the martingale 

£ • X t + 0„(w(t)) - M (w(O)) -n f ^(w(a)) da, 

■/ o 

and showing that for each fixed t, it is a Cauchy sequence in L 2 (Po) (and thus converges) 
as (j, tends to 0. This is achieved through spectral analysis, using the estimate (|2.10p . 
This estimate is obtained through a general argument of (anti-) symmetry (see the proof 
of [KV861 Theorem 4.1]), which has been systematized by [DFGW89] . Another difference 
is that the arguments of |KV86] are developped for general reversible Markov processes. 

3. Numerical approximation of the homogenized coefficients using the 

corrector equation 

3.1. General approach. Let £ G W 1 be a fixed unit vector. In order to approximate 
Ahom using the corrector (ft, we first replace the expectation in (|2.5p by a spatial average 
appealing to ergodicity: almost surely 

£-Aom£= hm / (£ + V^)-A(e + V0), (3.1) 

n ^°°Jq n 

where Q N = [0,N) d and f Q „ := N'^Q^- 

Recall that 4> is the solution to (|2.4p . which is a problem posed on Z rf . In order to turn 
(|3.ip into a practical formula, one needs a computable approximation </>jv of (j) on Qn, 
which would allow us to define an approximation An of ^hom by 

e • a n (, = / (e + v0jv) • ^(e + v^). 

Although ^hom is deterministic, this approximation An is a random variable. Provided 
the chosen approximation </>tv is "consistent", 

lim £ • ^at£ = £ • Aom£ 

almost surely. 

The starting point for a quantitative convergence analysis is the following identity: 

((£ • ^iv£ - t ■ A hom 2 ) = var [C ■ A N $ + (£ • (A N ) £ - £ • AomC) 2 . (3.2) 

The square root of the first term of the r. h. s. is called the random error or statistical 
error. It measures the fluctuations of An around its expectation. The square root of 
second term of the r. h. s. is the systematic error. It measures the fact that V^>at is only 
an approximation of V</> on Qn- In order to give quantitative estimates of these errors, 
we need to make assumptions on the statistics of A, and shall assume in the rest of this 
section that the entries of A have finite correlation length (that is, there exists Jz? c such 
that a(x, x + e,) and a(z, z + e^) are independent if \x — z\ ^ J*? c , for all i, j & {1, . . . , d}), 
which sligthly generalises the i.i.d. case and is enough for our purposes. 

In the following subsection, we introduce three approximation methods based on the 
above strategy, and recall the known (or expected) convergence rates for both errors. The 
third subsection is dedicated to a numerical study which completes the analysis in two 
respects: it allows to confirm/infirm numerically the expected convergence rates (if they 
are not precisely known), and it allows to make explicit the prefactors. 
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3.2. Methods and theoretical analysis. 

3.2.1. Homogeneous Dirichlet boundary conditions. The simplest approximation of the 
corrector <fi on Qn for N ^ 1 is given by the unique solution 4>n G L 2 (Qn H Z d ) to 

-V- A(e + V07v) =0 in Qjy, 

completed by the boundary conditions 4>n{x) = for all x G Z d \ Qtv- 

The approximation A at of Ah om is then defined by 



Z ■ A N C := / (f + V<^) • A(£ + Vc/w). 



As a direct consequence of homogenization, we have the almost sure convergence: 

Jim \A N - A hom | = 0. 

In terms of convergence rates, the starting point is identity (|3.2p . Although we do not 
have a complete proof, we expect the random error to scale as the central limit theorem, 
that is: 

var[e-^] 1/2 ~ N' d '\ (3.3) 
and the systematic error to scale as a surface effect (the corrector is perturbed on the 
boundary): 

\Z ■ (A N ) £ - £ ■ A hom £\ ~ N- 1 . (3.4) 
Let us give an informal argument for (j3.3|) in the case of ellipticity ratio (3/ a close to 1, 
that is, for A a perturbation of Id. Then, at first order, the equation for (fix takes the 
form 

-A<J) N = V • A£, 

completed by homogeneous Dirichlet boundary conditions on 8Qn- Let G denote the 
Green's function of the discrete Laplace equation on Qjv with homogeneous Dirichlet 
boundary conditions. Then, by the Green representation formula, 



4>n(x) = [ V y G(x,y)-A(y){;dy, (3.5) 
JQn 

V<Mx) = [ V x V y G(x,y)-A{y)tdy. (3.6) 
JQn 

We then appeal to the following spectral gap inequality (which is at the core [GOllj ): for 
any function A of a finite number of the i.i.d. random variables w e , we have: 

2\ 



var [X] ^ / ( sup 



UJ C 



dX 



doj f 



var [cj] , (3-7) 



where the supremum is taken w. r. t. the variable cj e , and var [uj] is the variance of the 
i.i.d. conductances. We shall apply this inequality to X = £ • Aat£. We first note that the 
weak formulation of the equation yields 

X = i-A N i = i £.A(£ + V^v), 
JQn 

so that by (|3.6p . we have for all e = (z, z + e{) 

IT- = -ftd&fc + V ^ N ^ + W I ^ ■ A(x)V x V Zi G(x, z)Cidx, 
e •> Qn 
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which, by symmetry of A and of G, we rewrite as 

I7 = Jfd&& + V iMz)) + j^£iV Zi J V x G(z,x) -A{x)idx. 
Using ()3.5p . this turns into 

dX 1 

It remains to take the supremum of this quantity w. r. t. uj e . In view of (|3.5|) . we have for 
all j 

oscVj(p(z) = sup Vj(f)(z) - inf V j<f>{z) ^ \V Xj V Vi G(x, z)\, 

which is bounded by a universal constant (in the discrete case the Green function is 
bounded). We thus have the desired variance estimate: 

var W < ^^<C(l + |V i ^(^)| 2 )>varM 

e 

J^C7(l + |Vi^(2r)| 2 )\var[w] 

= WwA dCNd + C ^ M ^Q^) vai[uj] 
< J_ 

since an elementary deterministic estimate yields 

|V^| 2 < P 2 \Qn\ = P 2 N d . 

The difficulty in the case of general ellipticity ratio is to treat the dependence of the 
Green's function with respect to u e , which yields additional terms of the form |V</>tv| 4 
which we do not control a priori (see |GNOa] in the case of periodic approximations), 
except for ellipticity ratio close to 1 (so that Meyers' estimate yields sufficient additional 
integrability of V</>jv). 

Note that the convergence rate of the random error (|3.3p is expected to depend on 
the dimension whereas the convergence rate of the systematic error (|3.4p does not. The 
combination of these (conjectured) estimates would yield the following convergence rate 
in any dimension 

3.2.2. Regularized corrector and filtering. In [GOllt fGQ12al IG112] . the following strategy 
was used. Instead of considering an approximation of the corrector <f>, we consider an 
approximation of the regularized corrector <j)^ of Lemma 12.51 The advantage of the reg- 
ularized corrector is that the Green's function associated with the operator fj, — V* • AV 
decays exponentially in terms of the distance measured in units of yU" 1 / 2 . Let N > 0. We 
denote by 4>^,n the unique solution in L 2 (Qn H 7L d ) to 

M0/i,iV - V • A(£ + Vc^tv) = in Q N , 

completed by the boundary conditions ^^{x) = for all x E r L d \Q^. 
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To define an approximation of Ah om , we introduce for all L £ N an averaging mask 
r)L '■ ^ d — > with support in Ql, and such that Jq l T]l = 1 and sup Z d |Vr/x| < L d+1 . 
For all iV ^ L > and /i > 0, we then set: 



J O.j 



(3.8) 



'Qz, 

We have the following almost sure convergence: 

lim lim \A^ NL -A hom \ 

Ai->0JV^L^oo 



0. 



Following the general approach, we have using the triangle inequality, and expanding 
the square: 

((£ • (A^, n ,l - A hom )0 2 ) 1/2 < <(£ • (\,n,l - A^ L )0 2 ) 1 ' 2 + var [A^ L ] 1 ' 2 + \A^ - A hom \, 

(3.9) 

where: 



C-A^lC := [ (e + V^)-^(e + V^)r/z 
JQt. 



(3.10) 



{•Apt -.= ((£ + V0 M ) • A(£ + V0 M )) . 

The first term is the error due to the fact that we replace </> M by the computable approx- 
imation 4>^n on Ql. As shown in |G112j the function (fiuN is a good approximation of 
c/)^ in the sense that we have the following deterministic estimate: there exists c > such 
that for all < L ^ N, we have almost surely 



N,L-A,, L \<^1- 



J\f\ d / 2 / AT \ d -V2 



L 



N 
N -L 



exp(-cVM(A r " L)). 



(3.11) 



This estimate essentially follows from the exponential decay of the Green's function mea- 
sured in units of /i" 1 / 2 . Hence, provided y/fl(N — L) 3> 1, the error between A^^^l and 
A^i is negligible. 

Within the assumption of finite correlation length, it is shown in |GQ11| that there 
exists q > (depending only on a and /3) such that 



var [£ • A^lQ 



d = 2 : L" 2 InV + ^ hiV, 
d > 2 : L~ d + /i 2 L 2 ~ d , 



(3.12) 



that is the central limit theorem scaling (up to a logarithmic correction in dimension 2) 
provided \xL < 1. (We expect this convergence rate to be optimal in general.) 
In |G()12a| (see also |GN()a| for d = 2), it is proved that 



d = 2 




d = 3 


/^ 3/2 , 


d = 4 


/i 2 1 In// 


d > 4 


/A 



(3.13) 
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Note that this scaling depends on the dimension and saturates at dimension d = 4. The 
proof of this estimate is interesting because similar arguments are used to analyze methods 
based on RWRE. It relies on an observation by one of us in [Moll] , and spectral calculus. 
Recall that for fixed unit vector £ £ M. d we denote by D = D* A£ the local drift and by 
the projection of the spectral measure onto fl. By definition, 

e-^e-e-^hom^ = ((e+D^-A^+D^g-^+D^-A^ + D^) 

= <(D^-D0).A(£ + D0) + (£ + D^)-A(D^-D0)). 

Using the weak form of the corrector equation (|2.8p . one has for all tp £ L 2 {Q), 

(D tp ■ A(£ + D 4>)) = 0. 

Taking tp = (j)^ — <j) v for v > 0, and using the symmetry of A, one obtains 

((D^-D^)-A(e + D0)) = = -<(£ + D0)-A(D<^-D<^)). 

Taking the limit v — > and using the strong convergence of D <p u to D (f> given by 
Lemma I2.7I this finally yields 

£-A^-£-A hom £ = ((D0 M -D0)- A(D0 M -D0)). 

Hence, by spectral calculus and taking the limit v — )■ in (|2.9|) . we obtain 

< i • A^ - A hom ^ = i? [ 1 & B (A). 

Jr+ 0" + A ) A 

This immediately implies that £ • A M £ — £ • Ah om £ ~ if A i— > A~ 3 is de a -integrable on M + . 
More precisely, what matters in this spectral integral is the behavior at the bottom of the 
spectrum. We cannot expect a spectral gap (which would hold in the periodic case since 
there is a Poincare's inequality on the torus) but we expect the bottom of the spectrum 
to be sufficiently thin to yield the result. It is convenient to rewrite the spectral integral 
as 

L (^W ea(A) * I! ¥TWx de » w + L \ de » ix) - 

By (|2.10|) the second term of the r. h. s. is bounded, and an elementary calculation shows 
that is a consequence of the following optimal estimate of the so-called "spectral 

exponents": for all < v ^ 1 and d ^ 2, 

/"^(A) < v d / 2+1 . (3.14) 
J o 

Suboptimal bounds with exponents d/2 — 2 were first obtained for d large in [Molll 
Theorem 2.4] using probabilistic arguments, optimal bounds up to dimension 4 (and up 
to some logarithm in dimension 2) in |GQ12aj using the spectral gap estimate (|3,7p and 
elliptic regularity theory, up to dimension 6 in [GM12aj by pushing forward the method of 
|GQ12aj (see |GM12bj for a generalization of this strategy which would yield the optimal 
exponents for all d > 2), and in any dimension in [GNOaj using the spectral gap estimate 
and parabolic regularity theory to obtain first bounds on the semi-group which turn after 
integration in time into bounds on the spectral exponents. 

Altogether, if we take N - L ~ N ~ L, and fi ~ iV -7 with 1 < 7 < 2 (so that (pTTTj) 
is of infinite order in N), the combination of the three estimates (|3.1ip . (|3.12p . and (|3.13|) 
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yields: 



\A/j,,n,L — ^4hom| 2 ) ^ 



d = 2 


jv- 


" 2 In 9 N, 


d = 3 


jv- 


-3 


d = A 


iV~ 




d > 4 


iV" 


-4 



(3.15) 



which yields the central limit theorem scaling up to dimension 4 at least. 

In dimension d > 2, since the fluctuations of A^l have the scaling of the central limit 
theorem (|3.12p . one may wonder whether the distribution of L d / 2 (£ • A^l^ — £ ■ (A^l) 
converges in law to a Gaussian random variable. Related results were obtained by Nolen 
[NoTT] . Rossignol [Rol2] . and Biskup, Salvi, and Wolff |BSW12| . 

3.2.3. Periodization method. The periodization method is a widely used method to ap- 
proximate homogenized coefficients, which consists in periodizing the random medium, 
see |Ow031 IEY07j . What is less known is that there may be several ways to periodize a 
random medium: 

• the periodization in space, 

• the periodization in law. 

Both periodization methods coincide in the specific case of i.i.d. conductances. The peri- 
odization in law is implicitly used in the very nice contribution [KFGMJj . 

Let A be a random conductivity function with finite correlation length. Let us begin 
with the periodization in space. It consists in approximating the corrector <f> on Qn by 
the QAr-periodic solution <p s ^ A with zero average to 

-V • A*> N (£ + V^ a ) = in Z d , 

where A& ,N is the Q jv-periodic extension of A\q n on that is for all k G 7L d and z G Qn, 
A&' N (kN + z) := A(z). The associated approximation A s ^ a of j4hom is then given by 



t ■ A s ^ :=i (£ + • A#> N {i + V# 

Jqn 



spa^ 



As a direct consequence of homogenization (see also |Ow03| ) , we have almost surely 

lim |^ a - A hom \ = 0. (3.16) 

In order to define the periodization in law, we have to make even more specific the 
structure of A. For all i G {l,...,d}, let gi : [a,/3] B — > [a,/3] be a measurable function 
depending only on a finite number of variables in B (recall that B is the set of edges), 
and let {u) e } ee B be a family of i.i.d. random variables in [a,/3]. We assume that the 
conductances {w e } eg B are given as follows: For all z G Z d and i G {1, . . . , d}, 

where t z O is the translation of O by z, i. e. for all e = (z',z' + e^) with z' G Z d and 
j G {1, ...,d}, T z cd e := oj( z+z i^ z+z i +e .y Since gi only depends on a finite number of 
variables, A has finite correlation length. The periodization in law consists in periodizing 
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the underlying i.i.d. random variables and then applying the deterministic function g. For 
all N ^ 1 we define oj*> n by: for all k G Z d , z G Q N , and i G {1, . . . , d}, u>t^ h+z+ei) := 
&(z,z+e t )> an d we se t for all z £ Z d 

A #jA r(z) := diag(g 1 (r z ui*' N ),...,g d (T z ui*' N )). 

Note that A# f jv is Q7v-periodic. It does coincide with A# ,N if gi(u) = Q(£d(Q }0+e .}) for 
some £?, in which case the conductances uj e are i.i.d. random variables, but not otherwise. 
We then consider the unique (^Ar-periodic solution cj}^ with zero average to 

-V-A #i7V (£ + V<^ w ) = OinZ d , 

and define 



e • 4re : = / a+ v^ w ) • a #i7V (£ + v<p 

JQn 

In order to prove the almost sure convergence 



N J 



lim |A| w -A hom | = 0, 

N-+00 



in view of (|3.16|) . it is enough to show that for large N, 

z ■ A s ^ - o(i) < e • a^^ c ■ a s ^c + o(i). 

By Meyers' estimates, there exists p > 2 such that for all N and almost surely, 

/ |V^T, f \V^r < N d . (3.17) 

JQn JQn 

Next we use the following alternative formula for A^ w (which holds by symmetry of A^jy): 

£ • A^ w £ = inf (/ (£ + V0) • A #)JV (£ + V0), is Q^v-periodic j . 

UQn ) 

Using (j) s ^ a as a test function and the fact that A^^r and A# ,N coincide on Qn-^ c , 
where «Sf c is the correlation-length of A, we get by Holder's inequality with exponents 
(p/2,p/(p-2)) 



< / (£ + vo-^(c+v^ a ) 

JQn 

= / (£ + v<^)^ # '^ + v<^ pa ) 

+ / (£ + V^ a )-(^#,iv-^ # ' JV )(e + V0 i 

< C-A # ^ + f3N- d [ (i + | V ^ a | 2 ) 

J Q m\Q m— r „ 



< £ • A^ + £iV- d |Qiv \ QA,-£ c | (p - 2)/p ||V^ a HiV) 
= £-A^ a £ + (l) 

using in addition (|3.17|) . The converse inequality is proved the same way. 
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Let us turn to convergence rates, and begin with the periodization in law, which is 
analyzed in |GNOa| IGNObj . The general approach takes the form: 



var 



As shown in |GNOal IGNObj , the random and systematic errors satisfy 



and 



var 



1/2 



< N -d/2 



(3.18) 



(3.19) 



The proof of (|3.19p is subtle. Let us give a hint of the proof in dimension d = 2, 3 for 
the case i.i.d. conductances. The idea consists in introducing a zero-order term as for the 
regularization, and we consider 

Z ■ tffrt := / (£ + V^) • A # , N (d + V^l), 
JQn 

where \i > 0, and </>Jf is the Qjv-periodic solution to 



V • A #)N (£ + V<^ w ) = 



m 



We also denote by A^ the approximation of A^ om defined in ()3.10p . Then by the triangle 
inequality, for all [i > 0, 



+ le • (4^) e - e • a^i + le • (^ - A om ) 



The first of the r. h. s. has the same scaling as the last term (repeating the arguments of 
spectral theory), that is f|3. 13|) . which is independent of N. The only term which relates 
N to ^ is the second term. In the i.i.d. case, we have by periodicity and stationarity 



Hence, since A l JT N {0) = A(0), 



Taw ^ 



f (e+v^)-A # ,^(e+v^ 

JQn 

(£ + V<^(0)) • 4>(0)(£ + V^(0)) 



(3.20) 



(e + v^(o)) • A(o)^ + v^(o)) - (e + v^(o)) • ^(o)(e + v^(o))) , 

and we only have to compare V0 M and at the origin. Using deterministic estimates 

on the Green's function, this yields 

If" £"^UI S ^exp(-c^iV), 

for some c depending only on a, /3, and d. Optimizing the error with respect to \x (taking 
for instance yf\x = ^^-L~ l ln(L/lnL)) yields the announced result (|3.19[) . 
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Let us make two comments on this proof. First, for d > 3 this strategy does not allow 
one to get (|3.19p since (|3.13p saturates at d = 4. Instead of comparing ^4hom to A^, we 
then compare ^4hom to a family of approximations A^^ for which we have 

£ • A k ^i - i ■ A hom £ < I — — de (A), 

for some p £ N which tends to infinity as k tends to infinity. This allows one to fully exploit 
(|3.14p . There are many possible choices for Ak tfl , as the one introduced in |GM12a] and 
defined by its spectral formula (and then translated back in physical space), and the one 
introduced in [GNOaj and defined by Richardson extrapolation in space. Second, the 
proof presented above crucially relies on the i.i.d. structure when writing ()3.20p . Indeed 
this relies on the fact that both statistical ensembles we use (the i.i.d. on B and on the 
set of edges of the A-torus) can be coupled (the ensemble on the torus being "included" 
in the ensemble on B). A similar construction can be made in the case of hidden i.i.d. 
variables as introduced above — but not in full generality. 

We expect these convergence rates to be optimal in general, so that the systematic 
error would scale as the square of the random error in any dimension (up to logarithmic 
corrections). As for the method with the zero-order term, we expect N d / 2 {^ ■ A N W £ — £ • 
(A l N w ) £) to converge in law to a Gaussian random variable, this time for all d ^ 2. 

In the case of the periodization in space, we still expect the random error to scale as 
the central limit theorem 

var[£.^] 1/2 < N~ d l\ (3.21) 
Yet, we do not expect the scaling of (|3.19p to hold in this case, and we rather conjecture 
that (unless the entries of A are i.i.d.) the systematic error scales as a surface effect (as 
for Dirichlet boundary conditions), namely 

\t-{AT)ti-ti- Ahamfl ~ N~\ (3.22) 

although we do not have a proof of this. The intuition behing this conjecture is that 
the imposed periodicity is not compatible with the underlying stationarity. A similar 
phenomenon occurs when considering a periodic problem and approximating the corrector 
on a domain which is not a multiple of the period and with periodic boundary conditions 
(so that the periodicity of the approximated corrector and that of the true corrector do 
not match), which yields an error which scales again as a surface effect (although the 
prefactor is usually much smaller than for Dirichlet boundary conditions). 

3.3. Numerical study. 

3.3.1. Homogeneous Dirichlet boundary conditions. We consider the simplest possible ex- 
ample: the case of Bernoulli variables. The conductances {w e } ee B are i.i.d. random vari- 
ables taking values a = 1 and /? = 4 with probability 1/2. 

In dimension d = 2, it is known that the homogenized matrix ^hom takes the form 
y/a/3Id = 2Id (the Dykhne formula, see [G112j for a proof). On Figure [1] we have plotted 

1 /2 

the estimates of the random error N i— > var [£ • A n£\ ' and of the systematic error A" i— >• 
|£ ■ (An) £ — £ • ^4hom£| in logarithmic scale. 

These errors are approximated by empirical averages of independent realizations, and 
intervals of confidence are given for the systematic error (corresponding to the empirical 
standard deviation) . The number of independent realizations in function of A^ is displayed 
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Figure 1. Dirichlet boundary conditions, d = 2, statistical error (red) 
rate 1.02 and prefactor 1.62, systematic error (blue) rate 0.98 and prefactor 
0.56. 



Table 1. 



N 


10 


20 


31 


56 


100 


208 


316 


Number of 
realizations 


1500 


6000 


14415 


47040 


150000 


648960 


1497840 



for completeness in Table [TJ As can be seen, the apparent convergence rates of the ran- 
dom error and of the systematic error are 1.02 and 0.98, respectively. This confirms the 
conjecture of (|3.3|) and (|3.4f) for d = 2. 

For numerical tests in dimension d = 3, we have to proceed slightly differently since 
there is no closed formula for the homogenized coefficient (the homogenized matrix is still 
a multiple of the identity by symmetry arguments). The approximation of the random 

1 /2 

error N i— > var[£-Ajv£] ' is unchanged. Yet, instead of plotting the systematic error 
N i — ^ |£ ■ (A N ) £ - £ ■ A hom £\, we plot an approximation of TV H> |f • (A N ) £ - £ ■ (A l ft w ) £| 
via empirical averages, where ^4^ w is the approximation of Ah om by periodization. In 
particular, in view of (|3. 19j) . we have by the triangle inequality 

|£-(A JV )£-£-^ hom £| > \C-(A N )Z-a.(A l ^)t\-CN- d ln d N 
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logic, (N) 

Figure 2. Dirichlet boundary conditions, d = 3, statistical error (red) 
rate 1.56 and prefactor 1.91, systematic error (blue) rate 0.90 and prefactor 
0.27. 

Table 2. 



N 


10 


15 


20 


25 


35 


43 


Number of 
realizations 


15000 


50625 


120000 


234375 


643125 


1192605 



for some C > 0, so that |£ • (A N ) £ - £ ■ A hom £\ and |£ • (A N ) £ - £ ■ (A l ft w ) f | are of the 
same order provided |£ • (An) ^ — £ ■ (yljy w )£| > N~ d ln d N (which we indeed observe 

numerically). These two errors N H> var • A N £\ 1/2 and N ^ |f • (Ajv) C ~ C • (-4jv w } C| 
are plotted on Figure [2] in logarithmic scale. The number of independent realizations in 
function of iV~ is displayed for completeness in Table As can be seen, the apparent 
convergence rate of the random error and of the (modified) systematic error are 1.56 and 
0.9, respectively. This confirms the conjecture of (|3.3|) and f)3.4j) for d = 3. 

3.3.2. Regularized corrector and filtering. We still consider the same simple two-dimensional 
example of the Bernoulli random variables taking values a = 1 and j3 = 4 with probability 
1/2, and for which the homogenized matrix is ^4hom = \fotf3I& = 2Id. In order to define 
^■ii,n,l completely, we need to choose L and (i in function of N, and define the averaging 
mask r]L- We have taken 
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' 1 1.5 2 2.5 

logio(N) 

FIGURE 3. Regularized corrector, d = 2, i .i. d. case, statistical error (red) 
rate 0.97 and prefactor 2.51, systematic error (blue) rate 1.55 and prefactor 
3.55. 



Table 3. 



N 


10 


17 


31 


56 


100 


177 


316 


Number of 
realizations 


1500 


4335 


14415 


47040 


150000 


510081 


1480338 



• L = 4N/5 and L = 3iV/5, 

• fj,= 125/iV 3 / 2 , 

• a piecewise affine mask, as plotted on Figure U] for N = 10 (the first one for 
L = 4N/5 and the second one for L = 3N/5). 

The theoretical predictions (|3. 12[) and ()3.13|) take the following forms with these parame- 
ters: 

var[e-^, L C] 1/2 < N^WN (3.23) 

1? • (A„,n,l) £ - £ • A hom £\ < N- 3 / 2 ln 2 N. (3.24) 

These two errors are plotted on Figure in logarithmic scale. The number of independent 
realizations in function of N is displayed for completeness in Table As can be seen, 
the apparent convergence rates for the random and systematic errors are 0.97 and 1.55, 
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FIGURE 4. Two masks for N = 10 

respectively. This shows the sharpness of the analysis. These numerical tests also give an 
idea on the prefactors in (|3.12p and (|3.13p . We observe that the systematic error decays 
faster than the random error, and that the prefactors are of the same order (roughly twice 
as big for the systematic error). 

The second series of tests on the regularization method aims at validating an approach 
which will be used to compare the two periodization methods. In particular, we consider 
now a two-dimensional case with correlations for which we do not have a closed formula of 
the homogenized coefficients. The statistics of the coefficients is defined as follows. We let 
{&( z ,z+ei)} z&fi ,ie{i,2} be i.i.d. variables following a uniform law in [0, 1]. We define W( 2j2 + ei ) 
to be a = 1 if for all z' such that \\z' — z\\oo ^ 2 we have u>r Z)Z + ei ) p (that is, if the 25 
hidden i.i.d. random variables are less than p), and Wf z>z + ei ) to be /3 = 4 if there exists z' 
with \\z' — zHoo ^ 2 such that u>r ZjZ+ei \ > p, where p is chosen so that p 25 = 1/2 (that is, a 
and (3 are equiprobable) . The typical realization of such conductances is made of islands 
of /3's of size 4 in a sea of a's. This example can be recast in the form of the correlated case 
described in Paragraph 13.2.31 Since the analysis of |GQ12al IGNOaj IGNObj also covers the 
case of finite correlation length, the theoretical predictions (|3.23p and (|3.24|) also apply to 
this case, provided we still take 

• L = 4N/5; 

• H = 125/iV 3 / 2 , 

• a piecewise affine mask, as plotted on Figured! 

Since we do not know j4hom a priori, we shall replace the systematic error by ./V i— >• 
|£ • (An t N t z — >1jv w ) CI j where A l ^ w is the approximation of j4hom by the periodization in 
law method. The combination of (|3.24p and (|3.19p indeed yields 

\C-(\,N,L-A^)^\ < iV- 3 / 2 ln 2 iV, 
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logio(N) 

Figure 5. Regularized corrector, d = 2, correlated case, statistical er- 
ror (red) rate 1 and prefactor 9.77, systematic error (blue) rate 1.35 and 
prefactor 20.9. 



Table 4. 



N 


10 


31 


100 


177 


246 


316 


Number of 
realizations 


1500 


14415 


150000 


469935 


907740 


1497840 



which we want to check numerically. This modified systematic error is plotted on Figure [5] 
in logarithmic scale. The number of independent realizations in function of N is reported 
on in Table HI As can be seen, the apparent convergence rate for the modified systematic 
error is close to 3/2, so that the true systematic error |£ • (A M) jv,l) — ^4hom£| has the same 
decay (up to possible logarithmic corrections) since 

ie • (a^ n>l ) - A hom z\ > le • (a^l - 4r) ei + cn- 2 in 2 n, 

for some C > due to (|3.19[) . Note that the asymptotic regime is more difficult to capture 
in the correlated case than in the i.i.d. since the typical lengthscale is 4 in the first case 
(the size of a typical island), and 1 in the second case. Yet these tests show it is possible 
to observe numerically a convergence with a rate larger than 1 in this correlated case. 
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logio(N) 

Figure 6. Periodization in law, d = 2, i.i.d., statistical error (red) rate 
1.01 and prefactor 1.48, systematic error (blue) rate 2.17 and prefactor 
1.29. 



Table 5. 



N 


10 


20 


31 


75 


135 


177 


Number of 
realizations 


1500 


14415 


498904 


469930 


780900 


1245782 



3.3.3. Periodization methods. In this subsection we first check numerically the sharpness 
of (|3.18p and (|3.19p on our simple two-dimensional Bernoulli case. On Figure [6] we have 

1 /2 

plotted the estimates of the random error N \- > var [£ • ^4^ w £] and of the systematic 
error N \— > |£ • (|yl^ w ) £ — £ • ^4hom£| m logarithmic scale. These errors are approximated 
by empirical averages of independent realizations, and intervals of confidence are given for 
the systematic error (corresponding to the empirical standard deviation). The number of 
independent realizations in function of N is displayed for completeness in Table [5j As 
can be seen, the apparent convergence rates of the random error and of the systematic 
error are 1.01 and 2.17 (note that the fluctuations are more important for the evaluation 
of the systematic error which is very small), which confirms the predictions (up to the 
logarithmic correction for the systematic error). In addition, the prefactors are again of 
the same order (slightly smaller for the systematic error), so that the systematic error is 
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Figure 7. Periodization in space, d = 2, correlated, statistical error (red) 
rate 0.99 and prefactor 6.17, modified systematic error (blue) rate 1 and 
prefactor 0.13. 



Table 6. 



N 


10 


31 


78 


100 


138 


177 


246 


316 


Number of 
realizations 


1500 


14415 


91260 


150000 


285660 


469935 


907740 


1497840 



really negligible in front of the random error. This will be further analyzed in the next 
subsection. 

The second series of tests deals with the correlated two-dimensional example introduced 
in Paragraph 13.3.21 In this case, we'd like to check numerically our conjecture (|3.22[) on 
the systematic error for the periodization in space method. As in the previous section, we 
replace the systematic error by the modified systematic error N \— > |£ • (A^ a — vl^ w ) £|. 
In view of (|3. 19j) , we indeed have 

ie • «^ a > - AomXi = le • (at - 4r) ^ + o(n- 2 in 2 n). 

1 /2 

On Figure [7] we have plotted the estimates of the random error N i— > var [£ • ^4jy w ^] and 
of the modified systematic error JV i->- |£ • (A^ w - A S T) CI in logarithmic scale. These 
errors are approximated by empirical averages of independent realizations, and intervals 
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of confidence are given for the systematic error (corresponding to the empirical standard 
deviation). The number of independent realizations in function of N is displayed for 
completeness in Table El As can be seen, the slopes of the random error and of the 
systematic error are 0.99 and 1, which confirms the conjectures (|3.2ip and (|3.22p on the 
periodization in space method. Yet it should be noticed that the prefactor of the systematic 
error is much smaller (around 40 times smaller) than the prefactor of the random error. 
This makes the systematic error hardly noticeable in practice. This would'nt be so clear 
in the case of dimension d ^ 3 (as on Figure [2]). 

It is instructive to compare the periodization in law A^^ of A to the periodization 
in space A N '& of A in terms of images. Typical realizations of these periodizations are 
represented on Figures [8] and [9] for N = 100 — 4 periods are reproduced. A close look at 
Figure [9] reveals the mismatch between stationarity and enforced periodicity (which does 
not appear on Figure [5] — this is the core of the "coupling" strategy) . 

3.3.4. Optimal numerical strategy. In this last paragraph, we propose a numerical strategy 
to obtain an approximation of the homogenized coefficients in the i.i.d. case at a given 
precision and at the lowest cost. In view of the previous subsections, it is clear that the 
periodization (in law) method minimizes both the random and systematic errors at N 
fixed. One feature we have not used yet is the scalability of the variance: if {^4j}i<ji<gfc are 
k independent realizations of a random variable A, then 

■ k 



var 



While empirical averages of independent realizations do not allow one to reduce the sys- 
tematic error, they do allow one to reduce the random error. This is particularly interesting 
since: 

• the dominant error is the random error, 

• it is computationally cheaper to solve several smaller linear problems than one 
single large problem (the solution cost of a linear system is always superlinear). 
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In view of the analysis of the periodization in law, in order to approximate homogenized 
coefficients within a tolerance 8 > 0, the cheapest computational way consists in solving 
kg independent periodic problems of size Ng, and take as approximation the empirical 
average of these kg realizations A^f. 



1 /ilaw 



(3.25) 



i=l 



where Ng is such that 
and kg such that 



8/2 = C syst N^ d ln d N 8 , 



8/2 



Crand-Ng 



-d/2 



where C sys t and C ranc i are the optimal prefactors in the estimates (|3.19|) and (|3. 18j) . re- 
spectively. Then we have 



/ 1 r 



2 V V2 



s: 8. 



To be precise, the associated computational cost is kgjd{Ng), where M h-> 7^(M) is the 
cost of solving a linear problem of order M (order of the matrix of the linear system) in 
dimension d. Since 7^ is superlinear, one readily convinces oneself that the best one can 
do is indeed (|3.25p . 

In the i.i.d. case of Figure El an approximation of these prefactors is 



C 



syst 



1.29, C rand = 1.48. 



4. Numerical approximation of the homogenized coefficients using the 

RWRE 

4.1. General approach. We now discuss how to approximate homogenized coefficients 
by simulating random walks. The simulation of a random walk has a very interesting 
feature: one does not need to generate a full environment a priori. Rather, it suffices to 
generate the environment along the trajectory of the random walk. This is particularly 
interesting in dimensions 3 and higher, where the random walk is transient, and visits only 
a vanishing fraction of the space. In fact, although the walk is recurrent in dimension 2, 
this last observation still holds in this case, since the time necessary to exit a box of size 
N is of order N 2 , while at such a time, the walk has typically visited a number of distinct 
sites of order iV 2 /ln(iV). 

In general terms, the strategy consists in simulating a large number of random walks, 
each in its own independent environment, and rely on Theorem 12.81 to recover the homog- 
enized coefficients. Keeping the environment fixed would be more difficult to analyse from 
a theoretical point of view, would certainly lead to larger errors (although we cannot prove 
this), and would force us to abandon the approach of generating the environment along 
the trajectory. 
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4.2. The discrete-time random walk. Although it is easier to see the link between the 
corrector equation and the continuous-time random walk, when it comes to simulations, 
it is more convenient to work with a discrete-time version of X, since there is no waiting 
times to compute. We define (Y n ) n& ^ to be the discrete-time Markov chain such that 
P^[X = z] = 1 and 



P" z [Yi = z'] 



p(z ~^ z') if z' ~ z, 

(4.1) 

otherwise, 



where p(z ~^ z') was introduced in f)2 . 13|) . Considering the "constructive" definition of the 
random variable X given in paragraph 12.3. H one can see that if Y n is the position of X 
after n steps, then (Y n ) n ^ is a discrete-time Markov process satisfying (|4.ip . As before, 
it will be convenient to consider the environment viewed by the walk Y, defined as 

w n = Qy n uj. 

Note that the probabilities p(z z') and p(z' z) need not be equal, hence it is not 
true in general that the counting measure is reversible for Y . A reversible measure ir for 
Y should satisfy, for any z,z' € 

ir(z) p(z z) = n(z ) p(z z). 

This relation holds if we choose tt(z) = Pu(z) (recall the definitions of Puj(z) and p(z zf) 
given in (|2.12p and ()2.13p respectively). This means that contrary to the continuous- 
time random walk, the discrete-time walk will preferably spend time on sites where Puj(z) 
is large. This effect has its counterpart concerning the environment viewed by Y: the 
measure IP is not invariant for (w n ) n 

£ N in general, but rather the following "tilted" version 

P, defined by 

d£(w) = dP(cA (4.2) 
(P) 

where for simplicity we write p(co) = Pui(0), and (p) = (p(ui)). In particular, the measure P 
has a density with respect to our initial measure P. Since p(oS) is never equal to or infinity, 
a property holds P-almost surely if and only if it holds with P-almost surely. 

We write Po for the measure PPq, and Eo for the associated expectation. Adapting 
slightly the arguments of paragraph 12. 3. 21 we obtain the following. 

Theorem 4.1 ( [KV86] ). Under the measure Po and as e tends to 0, the rescaled discrete- 
time random walk Y^ := (■^eY\ t / e \)tessL+ converges in distribution (for the Skorokhod 
topology) to a Brownian motion with covariance matrix 2^4^^, where 

^hom = (P)" 1 Aom = (2d K))" 1 A hom , (4.3) 

and ^4hom is a $ in (|2.5p . In other words, for any bounded continuous functional F on the 
space of cadlag functions, one has 



E 



F (y{>)) — -> E[F(B)], (4.4) 

e-s-0 



where B is a Brownian motion started at the origin and with covariance matrix lA^^. 
Moreover, for any ^ G l" 1 , one has 



(t-Y n y — 2£ ■ A<£ZZ- (4-5) 

n— >+oo 
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Before discussing numerical methods based on this theorem, we introduce some nota- 
tion. Let Y^ 1 ' , Y"( 2 ' , ... be independent discrete-time random walks evolving in the envi- 
ronments . respectively. We write Pg for their joint distribution, all random 



walks starting from 0, where tJ = (u!^\ uj^ 2 \ . . .). Tne environments , . . .) are 

themselves i.i.d. random variables distributed according to P, and we write P® for their 
joint distribution. We also write P® as a shorthand for the measure P®Pg. As usual, 
we simply replace "P" by "E" with the appropriate typography to denote corresponding 
expectations. 

4.3. Methods and theoretical analysis. In the two following paragraphs, we discuss 
numerical methods based respectively on the convergences in (j4.5|) and (|4.4p . We do so 
with the aim of explaining the quantitative estimates on the errors that can be proved, and 
the theoretical reasons lying behind the superiority of the method based on the convergence 
of the rescaled mean square displacement. 

4.3.1. Method based on the mean square displacement. As announced, we start with a 
method based on the convergence in (|4.5p . Recall that, by the definition of the tilted 
measure P in (14.21). we have 



The environments (u;W,u/ 2 ), 



n _1 E 



i 



n(p) 



E [p(w)(£-Y n ) 



By the law of large numbers, we know that for any fixed n, the quantity 

P (wW)(£ • Y n (1) ) 2 + • • • + P (w( fc ))(£ • y„ (fc) ) 2 



A,(n) 



kn {p) 



converges (almost surely) to the quantity in the right-hand side of equality 
to infinity. From the convergence in (|4.5p . we thus obtain 

OF ■ 4 disc t 



lim lim A^n) 

n— >+oo fc_ >+oo 



(4.6) 

(4.7) 
as k tends 

(4.8) 



The quantity Aj.(n) is what we will compute. It involves k random walks, each simulated 
in its own environment, up to time n. The expression involves the average {p) = 2d{u) e ), 
but this can be easily computed beforehand, so we assume that we have exact knowledge 
of this quantity. The convergence in (|4.8p can be made quantitative. 



Theorem 4.2 ! ( JM 121) ;. There exist constants q, C, c > such that for any k G N*, any 

e > and any n large enough, 

ke 2 ' 



A k (n) - 2£ • A^\ > (C/i d (n) + e)/n 



^ exp 



cn 1 



(4.9) 



where 



ln q n ifd = 2, 
1 i/d>3. 



This theorem ensures that it suffices to choose /c as a large constant times re 2 to guar- 



antee, with probability close to 1, that the difference between Ak{n) 
smaller than C/n (or C In 9 (n)/n if d = 2) for some constant C. 



and 2£ • Af^t. is 
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We decompose the error 

\A k (n) - 2£ ■ A*ZZ\ 

into the sum of 

A k {n)-n- l t Q \(£-Y n ) 2 



+ 



(4.10) 



We call the first difference the statistical error, and the second difference the systematic 
error. The theorem is proved if we show the following two inequalities: 



AfcCnJ-n^to (£ • Y n ) 2 > e/n 



he 2 

^ exp ( - — 



err 



n _1 En 



4 disc a 



Cji d {n)/n. 



(4.11) 
(4.12) 



Inequality (|4.1ip is controlled using classical large deviations theory, noting that A k {n) is 
a sum of i.i.d. random variables. In the same vein, it is also possible to show (see |GM12b[ 
Proposition 5.1]) that for any sequence k n tending to infinity with n, 



k n I A kn (n 



n 



(distr.) 



>A/"(0,v) 



(4.13) 



where M(0, v) is a Gaussian random variable of variance v, and 



i 



m-A^o 2 - 



(4.14) 



Although the left-hand side of inequality (|4.12p is deterministic, its proof requires a more 
subtle analysis. The problem can be rephrased by saying that a quantitative control of 
the convergence in (|4.5[) is needed. 

For simplicity, we will instead discuss how to obtain a quantitative control of the con- 
vergence in (|2.20p . that is, for the continuous-time random walk. Recall that in para- 
graph [2321 w e decomposed £ -Xt into Mt + Rt, where Mf is a stationary martingale under 
the measure ¥q, and Rt is a remainder. The martingale property and the stationarity of 
the increments guarantee that Eo[M t 2 ] grows linearly with t, and in fact (see (|2.24|) and 

E [M 2 ] = 2t £ ■ A hom £. (4.15) 

We stress that the fact that there is an equality above, and not just an asymptotic equiva- 
lence, is in itself a "miracle" : it says that for any t, the expectatation of M 2 is exactly the 
expectation of (£ • Bt), where B is the limiting Brownian motion (and moreover, note that 
we did not need any quantitative information whatsoever when deriving this identity). 
Using (|4,15p and the fact that £ • Xt = Mt + Rt, we thus obtain 

E [(£ • X t ) 2 ] = 2tC ■ A hom Z + E [R 2 ] + 2E [M t R t ] ■ (4.16) 

In paragraph 12.3.21 we have sketched the argument leading to a proof of the fact that 
Rt/Vi tends to in L 2 (Fq) (see (|2.27p ): the problem is reduced to a proof of the fact 
that the spectral measure is light enough in the neighborhood of to integrate A -1 (see 
(I2.1(jp ) ; and this in turn is proved via a general argument of symmetry. 

This general argument is interesting, and covers arbitrary stationary distributions of 
conductances. However, it gives no information on the rate at which i -1 Eo [R 2 ] converges 
to 0, and this is inherent to such a high level of generality. When the conductances 
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are independent (or typically with finite correlation- length) , we recall that the following 
optimal estimates of (|3.14|) hold (see |GNOaj ): for all v > 0, 



f de a (A)<z// 2+1 . 
J o 



In view of (|2.27|) . such a control of the behavior of the spectral measure at the edge of the 
spectrum should give information on the speed of convergence of t-% \ Rf\ to 0. Indeed, 
one can show using basic calculus manipulations (see the proof of [Molll Proposition 8.2]) 
that (|3.14p implies that in fact, Eo [Rt] remains bounded as t tends to infinity (or grows 
no faster than ln(i) in dimension 2). The crux of the proof is thus to show that (|3.14p . 
In the case of discrete time, our estimates of the spectral are slightly weaker (we have 
a logarithmic divergence in dimension 2, and optimal estimates up to dimension 6), see 
|GM12al Appendix A]. This is however sufficient to prove Theorem 14.21 

In the decomposition (|4.16p . there finally remains to study the cross-product Eo [MtRt]. 
Since we know that roughly speaking, Eo [Rt] remains bounded as t tends to infinity, a 
naive Cauchy-Schwarz estimate would ensure us that |Eo [MtRt] I grows no faster than yt. 
But in fact, a second "miracle" occurs here, since this cross product turns out to be ! 
Recall that it can be written as 

E [(£ • X t + 4>{X t ,u) - 0(0, oj)) (<f>(X t ,uj) - 0(0, w))] . (4.17) 

To see that this expectation is 0, it suffices to show that its derivative at time vanishes, 
thanks to the Markov property and to the stationarity of (u(t))t^o- Using (|2.11|) . we 
obtain that this derivative is equal to 

/X> 0j2 (£ • z + 0(z,oj) - 0(O,u/)) (<}>(z,gj) - <t>(0,u))\ = {A(£ + V0) • V0> , 

and the latter is equal to by (j2.4[) (the fact that (|4.17|) vanishes is also clear from the 
alternative construction of the corrector based on orthogonal projections, as was done for 
instance in |MP07j ). 

This completes our sketch of the proof of the continuous-time analog of Theorem 14.21 
For discrete time, we refer the reader to [GM12bj . 

4.3.2. Methods based on other functions of the random walk. We now turn our attention to 
numerical methods based on the convergence in (|4.4p . and in fact, we continue to discuss 
instead the continuous-time analog (|2.19p . We focus on a special case of (|2,19p . namely, 
the fact that for any bounded continuous / : R d — > R, 



Er 



/ 



(Xt 

Wt 



f— >+oo 



> E[f(Bi)], (4.18) 



where B\ is a Gaussian vector with covariance matrix 2^4hom- In the preceeding chapter, 
we highlighted two "miracles" that made the difference between the l.h.s. of (|2.20p and its 
limit as small as t -1 . Both these miracles are very special of the square function. Although 
we cannot prove it, for a generic functional, we expect the error to be of order at least 
t" 1 / 2 in general. We now present the known upper bounds on the error. 



Theorem 4.3. Let B\ be a Gaussian vector with covariance matrix 2A\ lom . 
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(1) [Mol2c] /// is smooth and bounded, then for any 5 > 0, there exist q,C such that 

t- 1 ^ ifd=l, 



En 



/ 



Wt 



E[f(B 1 



^ C 



ln q (t) t' 1 ^ 
£-1/2+6 

(2) [Mol2bj There exists q,C such that for any f G R d , 



ifd = 2, 
3. 



(4.19) 



sup 



Vt 



€ x 



Vt 



^ X 



< c 



t -i/w 


i/d = l, 


in^t) r x /io 


i/d = 2, 




ifd = 3, 


r l/5 


ifd^A. 



(4.20) 



Part (2) of the theorem corresponds to considering the function /(z) = lc. z ^. x , and is 
often known as a Berry-Esseen estimate. 

For comparison, note that if Xt was replaced by a sum of i.i.d. centred random variables 
with finite third moment and covariance matrix 2A} lom , then the left-hand sides of (|4.19p 
and (|4.20p would both be 0(t~ l l 2 ) (see for instance |Pet Theorem 5.5]). It was recently 
realized (but yet to be written) that the bound in (|4.19p for d = 2 can be improved to 
t~ l / 2+5 for any 6 > 0. The conjectured optimal rates in this inequality are 



t"V4 if d = 1, 

ln(i) t~ 1 / 2 if d = 2, 
if O 3. 



(4.21) 



Remark 4.4. To be precise, if Xt was replaced by a sum of i.i.d. centred random variables, 
say Z\,...,Zt (t E N here) , then the characteristic function of the normalized sum would 
satisfy (for any ( G K) 



E 



exp i£ 



E 



exp^jj =exp(-^-,i 7r l + 0(^)j, 



where we assumed that Z\ has finite moments of order 4. The first term in the expansion 
gives the characteristic function of the limiting Gaussian random variable, and there is 
indeed a correction of order i -1 / 2 . But this correction vanishes in the special case when 
KZf = 0, and in particular if the distribution of Z\ is invariant under the transformation 
z I—)- —z. Numerical simulations suggest that similar cancellations sometimes also happen 
for the random walk in random environment, as we will discuss below. 



Without entering into details, we now sketch the main steps of the proof of Theorem[ 
Given that, roughly speaking, En[-R 2 ] remains bounded, it is not too hard to control the 
error due to the absence of the "miracle of orthogonality" . We thus turn our attention to 



CONVERGENCE RATES IN STOCHASTIC HOMOGENIZATION 



35 



the absence of the first mentionned "miracle" , which consists in the fact that the relation 
(|4.15|) is exact, regardless of whether the martingale Mt is actually closely resembling a 
Brownian motion or not. For a generic function /, it is no longer true that the difference 



E 



J \Vt 



E[f(Bi)] 



(4.22) 



vanishes, and the proof of Theorem 14.31 requires that we find bounds on this quantity. 
Recall that we derived the fact that the difference in (|4.22p converges to from the 
convergence of t~ 1 Vt in (|2.26p . The two main steps of the proof of Theorem 14.31 are: 

(i) to estimate the speed at which the variance of t~ l Vt tends to (i.e. to get an L 2 
estimate on the speed of convergence in (|2.26p ). and 

(ii) to rely on a general quantitative version of the central limit theorem for martingales 
which transports the estimate obtained in (i) into an upper bound on (|4.22p . 

At least for d ^ 4, our control of the variance of t~ l Vt is optimal, and the reason for the 
not-so-intuitive exponent 1 /5 appearing in part (2) of Theorem 14.31 is hidden in part (ii) 
of the proof, that is, in the general quantitative central limit theorem for martingales. 
Surprisingly, this general result is however optimal |Mol2a] , Yet, we conjecture that the 
exponents obtained in part (2) of Theorem l4.3l are not optimal, and that the rates obtained 
in part (1) of the theorem may in fact hold as well in part (2). 



4.4. Numerical study. For our numerical tests, we continue with the case introduced in 
paragraph 13.3.11 of i.i.d. conductances (cj e ) ee B with 

P [u e = 1] = P [w e = 4] = 1/2. 

Recall that by symmetry considerations, we know that in this case, ^ s a multiple of 

the identity. For the two-dimensional case, we know moreover that ^4hom = 2 Id, and thus 
by ()4.3p . that 4^ = 1/5 Id. We do not know of such a formula when d = 3, and for the 
purpose of the analysis of the methods, we may use an approximation of the homogenized 
coefficients obtained with the periodization method of paragraph 13.2.31 

For the simulation of the random walk, we generate the environment only along the 
trajectory of the walk. More precisely, when the walk arrives at some point, it first checks 
if the neighboring conductances have already been "seen". If so, then the value already 
generated is used, else a value is taken at random and stored for future usage. Roughly 
speaking, the walk up to time n discovers of order n conductances, independently of the 
dimension d ^ 2 (except for a logarithmic correction in dimension 2). 

The methods based on simulating random walks have two main interesting features. 
One is that their effectiveness is fairly insensitive to the dimension. The other is that 
the computation is massively parallel, since it requires to simulate a lot of independent 
random walks, each in its own independent environment. 



4.4.1. Method based on the mean square displacement. We start by investigating numeri- 
cally the method based on the computation of the mean square displacement of the random 
walk. As for the methods based on the corrector, we investigate separately the statistical 
and systematic errors, that is, the two terms in the sum (|4.10p . For the statistical error, 
we focus on 

2" 



(i fc (?i)-n- 1 E [(£-Y n ) 



Var 



A k {n) 
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n 


10 


20 


40 


80 


160 


320 


640 


1280 


oo 


Variance 


0.40 


0.39 


0.380 


0.373 


0.369 


0.367 


0.3653 


0.3647 


0.3632 



Table 7. Numerical estimates for (|4.24p and the theoretical limiting value 
in dimension 2. 



n 


10 


20 


40 


80 


160 


320 


640 


1280 


oo 


Variance 


0.20 


0.19 


0.188 


0.186 


0.184 


0.184 


0.1836 


0.1835 


0.1829 



Table 8 . Numerical estimates for (|4.24|) and the theoretical limiting value 
(computed with a numerical estimate for ijj") in dimension 3. 



where we write Var [•] for the variance with respect to the measure 
sum of independent random variables, we have 

2" 



Var 



A k (n) 



1. 

k 



-Var 



pM (t-Y r 

(p) 



Since Ak(n) is a 



(4.23) 



A simple variation of the proof of the convergence in (|4.13p shows that 

lim Var 



n— >+oo 



pM (t-Yn 



(p) 



where v was defined in (|4.14p . Table [7] shows our numerical estimations in dimension 2 for 



Var 



(P) 



(4.24) 



for several values of n, and compares it to the predicted limiting value. As we see on 
(|4.23|) . the variance of our estimator Ak(n) is then obtained by dividing this value by the 
number k of walks we run. Table [5] presents the same results for the 3-dimensional case. 

Turning to the systematic error, the theoretical prediction is that displayed in (|4.12|) . To 
test this, we computed Ak n (n) with k n = K{n)n 2 , where K{n) is some large number. This 
choice of k n ensures that the random fluctuations are washed out, so that Ak n {n) is very 
close to its expectation. In practice, we chose K{n) = 10 4 for n 320, and K{n) = 10 3 
for larger values. 

For the two-dimensional case, our results are reported on Figure [TUJ In this case, the 
theoretical prediction in (|4.12p contains a polylogarithmic correction to the rate of decay 
in 1/n. We conjecture that the true rate is actually ln(n)/n. On the log-log plot, this 
would create a correction to the expected slope of —1 of the order of l/ln(n), which is 
between 0.17 and 0.14 for n between 320 and 1280, roughly corresponding to the observed 
slope. 

In three dimensions, our findings are reported on Figure [TTJ Here, the observations 
match very well with our predicted convergence rate of 1. 
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4.4.2. Methods based on other functions of the random walk. We now turn our attention 
to methods based on computing other functions of the final position of the random walk 
than the squared displacement. For any reasonable function /, we can devise an estimator 
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as we did for the special function f(x) = (£ • x) 2 , which takes the form 



A' 



^ W := W) E^ (4) )/( y n W /V^)- (4-25) 



We have 

if 



E®[A{(n)]=E [f(Y n /^% 
which (for suitable /) converges to E[f(Bi)], where B% is a Gaussian vector with covariance 
matrix 2-4??°. Since this more general estimator is still a sum of i.i.d. random variables, 
the statistical error is easy to understand, as we have 



Var 



AUn) 



^Var 
k 



(p) V V™ 



where the last variance converges (for suitable /) to a constant that can be explicitely 
written. We now focus our numerical investigation to the systematic error, that is, the 
difference between Ko[f(Y n /y/n)] and the limiting value E[f{B\)\. Based on the analogy 
with the continuous-time case, we expect that the systematic error will decay as n~ 1 / 2 for 
d ^ 2, possibly with a logarithmic correction in dimension 2. 
In order to test this prediction, we chose the function 

f{x) = exp (-^p) ■ (4.26) 

In the case we investigate numerically, the conductances being i.i.d., we know that the 
homogenized matrix is of product form, say = a 2 Id, and actually we also know 

that a 2 = 2/5. A simple computation shows that in this case and for / defined in (|4.26p . 
the limiting value is given by 

E[f(B 1 )\ = (a 2 + l)^ 2 , 

which is equal to 5/7 in dimension 2, and which we can approximate thanks to our previous 
numerical estimation of a 2 in dimension 3. 

The systematic errors we obtained numerically are shown on Figures [12] and [13] for 
dimensions 2 and 3 respectively. 

Surprisingly, the observed convergence rates are far better than the predicted ones, 
being close to 1 in both cases instead of the predicted 1/2. Similar rates where observed 
for other choices of the function /. 

We believe that these surprising rates are due to special cancellations which do not hold 
for arbitrary distributions. As was explained already in Remark 14.41 such cancellations 
happen already in certain cases when one considers sums of i.i.d. random variables, in 
particular when the distribution of the random variables is invariant under the transfor- 
mation z i — y —z. Another way to guess what the correct order of the correction is to 
compute an asymptotic expansion in terms of the ellipticity ratio, in the spirit of what 
was done in [GQ11[ Appendix]. If one looks at the first order correction, it turns out 
that it behaves as t^ 1 ! 2 in general, but cancellations occur if the distribution is symmetric 
under the transformation z i— > —z, and the correction becomes of order t . In view of 
these two facts, we believe that the convergence rates obesrved on Figures [P21 and [TBI are 
somehow miraculous, and due to the symmetry of the distribution of the conductances 
under the transformation z i— >■ — z. 
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logio(N) 

Figure 12. The systematic error in two dimensions for / as in (|4.26p . rate 
1.09 and prefactor 0.14. 




logio(N) 

Figure 13. The systematic error in three dimensions for / as in (|4.26p . 
rate 1.06 and prefactor 0.07. 
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-O Edge with conductance equal to 1 



Periodic cell 



Figure 14. The periodic environment, without z i— > —z symmetry 
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Figure 15. The systematic error in the periodic environment for f(x, y) 
sin(x), rate 0.49 and prefactor 0.12. 



In order to convince ourselves that the convergence rates observed are indeed non- 
typical, we reverted to a study of a periodic environment, choosing it so that it is not 
symmetric under the transformation z i— > —z. A periodic environment should provide bet- 
ter convergence rates than any generic and truly random environment (in any dimension) . 
We chose a simple 3-periodic cell displayed on Figure HH We chose f(x,y) = sin(x), so 
that Ko[f(Y n /y/n)] tends to as n tends to infinity (and we thus do not need to compute 
the homogenized matrix to study the systematic error). The results we obtained are shown 
on Figure H5J As expected, they display a convergence rate of 1/2. 
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Figure 16. The systematic error in two dimensions for / as in (|4.27p and 
z = 1/2 (red), rate 0.8 and prefactor 0.2, or z = 1/4 (blue), rate 1.0 and 
prefactor 0.7. 



Finally, we investigated the convergence rates for the non-smooth function 

f(x) = l f . I<2) (4.27) 

where we took £ to be the first vector of the canonical basis and z = 1/2 or z = 1/4. For 
the systematic error, Theorem 14.31 ensures a convergence rate of 1/10 in dimension 2, and 
of 1/5 in dimension 3, up to logarithmic corrections. We believe that the exponent can be 
pushed to 1/5 in dimension 2 with a refined argument (yet to be written), but the proof 
cannot be pushed to higher exponents |Mol2aj . 

The main goal of these numerical investigations is to see whether this limitation to the 
exponent 1/5 is an artefact of the method of proof. Our results are displayed on Figure [TBI 
for d = 2, and on Figure [T71 for d = 3. 

While the onset of a nice asymptotic regime is delayed, the results suggest that indeed 
the convergence rates ultimately settle to a behavior similar to that observed on Figure [T2l 
or 113} that is, ultimately displaying a convergence rate close to the value 1. 



5. CONCLUSION 

In this article we have recalled qualitative and quantitative results on stochastic homog- 
enization of discrete linear elliptic PDEs and of random walks in random environments 
on 7j d . This has allowed us to make a rather complete picture of numerical methods to 
approximate homogenized coefficients, based both on the corrector equation and on the 
random walk. Numerical tests have confirmed the sharpness of the analysis, supported 
some conjectures, and put some interesting phenomena into evidence (such as ungener- 
ically fast decay rates due to specific symmetries of the environement). We hope this 
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Figure 17. The systematic error in three dimensions for / as in (|4.27|) 
and z = 1/2 (red), rate 0.9 and prefactor 0.2, or z = 1/4 (blue), rate 1.0 
and prefactor 1.0. 



contribution will help mathematicians identify challenging conjectures, and help practi- 
tioners make mathematically-based choices on the method to use in more concrete cases. 

We have only considered discrete elliptic equations. The discrete case has the advantage 
of being numerically inexpensive to simulate compared to the case of continuum linear 
elliptic equations (for which we have to appeal to approximation methods such as the 
finite element, finite difference or finite volume methods). The extension of the results 
of this paper to the continuum case is currently under investigation, and |GQ1H IG012a| 
have already been extented to the continuum case |GQ12bj . The extension of the RWRE 
approach to the continuum space is yet to be done. We hope to address these issues in 
future works. 
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